Máquina de bolas - análisis

Resumen Ejecutivo

Este análisis exploratorio examina datos operacionales de molinos de bolas del caso propuesto, enfocándose en la identificación de patrones para mantenimiento predictivo y optimización operacional. Los datos contienen 130k registros con 64 variables que abarcan desde parámetros operacionales hasta indicadores de fallas.

Objetivos principales:

1. Identificar patrones precursores de fallas en equipos críticos
2. Analizar relaciones entre variables operacionales y eficiencia
3. Proporcionar insights para optimización del proceso de molienda

1. Configuración y Carga de Datos

# Importación de librerías

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import warnings
from scipy import stats
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
import datetime as dt

from utils import molinos_data

# Configuración de visualización
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")
warnings.filterwarnings('ignore')

# Configuración de plotly
import plotly.io as pio
pio.templates.default = "plotly_white"

print("✅ Librerías cargadas exitosamente")

# Carga de datos
df = molinos_data()

# Conversión de tipos
df['timestamp'] = pd.to_datetime(df['timestamp'])
df['falla_en_7d'] = df['falla_en_7d'].astype(bool)
df['falla_en_14d'] = df['falla_en_14d'].astype(bool)
df['falla_en_30d'] = df['falla_en_30d'].astype(bool)

print("✅ Datos cargados y procesados")
print(f"📊 Shape del dataset: {df.shape}")
print(f"📅 Rango temporal: {df['timestamp'].min()} a {df['timestamp'].max()}")

✅ Librerías cargadas exitosamente
✅ Datos cargados y procesados
📊 Shape del dataset: (131478, 64)
📅 Rango temporal: 2023-01-01 00:00:00 a 2025-07-02 00:00:00

Información General del Dataset

# Información general
print("=== INFORMACIÓN GENERAL DEL DATASET ===")
print(f"Registros: {len(df):,}")
print(f"Variables: {len(df.columns)}")
print(f"Memoria utilizada: {df.memory_usage(deep=True).sum() / 1024**2:.2f} MB")

# Distribución por molino
print("\n=== DISTRIBUCIÓN POR MOLINO ===")
molino_dist = df['molino_id'].value_counts().sort_index()
print(molino_dist)

# Distribución temporal
print("\n=== DISTRIBUCIÓN TEMPORAL ===")
df['year_month'] = df['timestamp'].dt.to_period('M')
temporal_dist = df['year_month'].value_counts().sort_index()
print(f"Período más antiguo: {temporal_dist.index.min()}")
print(f"Período más reciente: {temporal_dist.index.max()}")

# Valores faltantes
print("\n=== VALORES FALTANTES ===")
missing_data = df.isnull().sum()
missing_pct = (missing_data / len(df)) * 100
missing_summary = pd.DataFrame({
    'Valores_Faltantes': missing_data,
    'Porcentaje': missing_pct
}).sort_values('Porcentaje', ascending=False)

if missing_summary['Valores_Faltantes'].sum() > 0:
    print(missing_summary[missing_summary['Valores_Faltantes'] > 0].head(10))
else:
    print("✅ No se encontraron valores faltantes")
# Valores
print("\n=== VARIABLES ===")
print(df.info())

=== INFORMACIÓN GENERAL DEL DATASET ===
Registros: 131,478
Variables: 64
Memoria utilizada: 81.29 MB

=== DISTRIBUCIÓN POR MOLINO ===
molino_id
M1    21913
M2    21913
M3    21913
M4    21913
M5    21913
M6    21913
Name: count, dtype: int64

=== DISTRIBUCIÓN TEMPORAL ===
Período más antiguo: 2023-01
Período más reciente: 2025-07

=== VALORES FALTANTES ===
                      Valores_Faltantes  Porcentaje
temperatura_trend_7d                138    0.104961
vibracion_trend_7d                  138    0.104961
throughput_trend_24h                 66    0.050199
energia_trend_24h                    66    0.050199

=== VARIABLES ===
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 131478 entries, 0 to 131477
Data columns (total 65 columns):
 #   Column                          Non-Null Count   Dtype         
---  ------                          --------------   -----         
 0   timestamp                       131478 non-null  datetime64[ns]
 1   molino_id                       131478 non-null  object        
 2   turno                           131478 non-null  object        
 3   feed_rate                       131478 non-null  float64       
 4   velocidad_rotacion              131478 non-null  float64       
 5   velocidad_porcentaje_critica    131478 non-null  float64       
 6   nivel_carga_bolas               131478 non-null  float64       
 7   densidad_pulpa                  131478 non-null  float64       
 8   agua_adicionada                 131478 non-null  float64       
 9   presion_ciclones                131478 non-null  float64       
 10  vibracion_cojinete_feed_h       131478 non-null  float64       
 11  vibracion_cojinete_feed_v       131478 non-null  float64       
 12  vibracion_cojinete_discharge_h  131478 non-null  float64       
 13  vibracion_cojinete_discharge_v  131478 non-null  float64       
 14  vibracion_shell_h               131478 non-null  float64       
 15  vibracion_shell_v               131478 non-null  float64       
 16  vibracion_pinion                131478 non-null  float64       
 17  vibracion_gearbox               131478 non-null  float64       
 18  temp_cojinete_feed              131478 non-null  float64       
 19  temp_cojinete_discharge         131478 non-null  float64       
 20  temp_aceite_lubricacion         131478 non-null  float64       
 21  temp_motor_principal            131478 non-null  float64       
 22  temp_gearbox                    131478 non-null  float64       
 23  corriente_motor                 131478 non-null  float64       
 24  potencia_activa                 131478 non-null  float64       
 25  voltaje_motor                   131478 non-null  float64       
 26  factor_potencia                 131478 non-null  float64       
 27  presion_aceite_principal        131478 non-null  float64       
 28  flujo_aceite                    131478 non-null  float64       
 29  nivel_tanque_aceite             131478 non-null  float64       
 30  calidad_aceite_ppm              131478 non-null  float64       
 31  consumo_energetico_especifico   131478 non-null  float64       
 32  throughput_real                 131478 non-null  float64       
 33  eficiencia_molienda             131478 non-null  float64       
 34  granulometria_producto_p80      131478 non-null  float64       
 35  nivel_desgaste_liners           131478 non-null  float64       
 36  horas_operacion_acumuladas      131478 non-null  int64         
 37  ciclos_arranque_parada          131478 non-null  int64         
 38  carga_circulante                131478 non-null  float64       
 39  eficiencia_clasificacion        131478 non-null  float64       
 40  work_index_bond                 131478 non-null  float64       
 41  dureza_mineral                  131478 non-null  float64       
 42  humedad_mineral                 131478 non-null  float64       
 43  granulometria_feed_p80          131478 non-null  float64       
 44  densidad_mineral                131478 non-null  float64       
 45  contenido_arcillas              131478 non-null  float64       
 46  abrasividad_ai                  131478 non-null  float64       
 47  temperatura_ambiente            131478 non-null  float64       
 48  humedad_relativa                131478 non-null  float64       
 49  falla_en_7d                     131478 non-null  bool          
 50  falla_en_14d                    131478 non-null  bool          
 51  falla_en_30d                    131478 non-null  bool          
 52  tipo_falla                      131478 non-null  object        
 53  severidad_falla                 131478 non-null  int64         
 54  dias_hasta_falla                131478 non-null  float64       
 55  vibracion_trend_7d              131340 non-null  float64       
 56  temperatura_trend_7d            131340 non-null  float64       
 57  energia_trend_24h               131412 non-null  float64       
 58  throughput_trend_24h            131412 non-null  float64       
 59  ratio_p80_feed_producto         131478 non-null  float64       
 60  potencia_especifica_neta        131478 non-null  float64       
 61  eficiencia_energetica_teorica   131478 non-null  float64       
 62  anomaly_score_vibration         131478 non-null  float64       
 63  anomaly_score_electrical        131478 non-null  float64       
 64  year_month                      131478 non-null  period[M]     
dtypes: bool(3), datetime64[ns](1), float64(54), int64(3), object(3), period[M](1)
memory usage: 62.6+ MB
None

# Información general
df.head(10)

	timestamp	molino_id	turno	feed_rate	velocidad_rotacion	velocidad_porcentaje_critica	nivel_carga_bolas	densidad_pulpa	agua_adicionada	presion_ciclones	...	vibracion_trend_7d	temperatura_trend_7d	energia_trend_24h	throughput_trend_24h	ratio_p80_feed_producto	potencia_especifica_neta	eficiencia_energetica_teorica	anomaly_score_vibration	anomaly_score_electrical	year_month
0	2023-01-01 00:00:00	M1	A	255.526011	13.877000	76.937236	32.085144	75.608893	82.003212	98.694143	...	NaN	NaN	NaN	NaN	86.723185	16.267886	0.758328	0.431301	0.961619	2023-01
1	2023-01-01 00:00:00	M2	A	293.793464	13.789087	76.449823	32.027222	74.903245	88.460830	89.432729	...	NaN	NaN	NaN	NaN	104.352346	15.111197	0.872687	0.268677	0.108372	2023-01
2	2023-01-01 00:00:00	M3	A	296.830860	13.874538	76.923581	33.838441	74.000979	86.905074	122.120387	...	NaN	NaN	NaN	NaN	96.711104	18.555213	0.684239	0.185597	1.047968	2023-01
3	2023-01-01 00:00:00	M4	A	268.552994	13.450880	74.574725	29.904668	73.309069	59.344657	109.740612	...	NaN	NaN	NaN	NaN	104.111652	17.422468	0.752651	0.231585	0.053925	2023-01
4	2023-01-01 00:00:00	M5	A	230.674766	13.956928	77.380372	32.439938	79.171183	69.664273	77.401938	...	NaN	NaN	NaN	NaN	96.249339	19.060690	0.695311	0.584719	0.905303	2023-01
5	2023-01-01 00:00:00	M6	A	279.021934	13.624530	75.537482	33.065775	69.762263	81.290097	109.618477	...	NaN	NaN	NaN	NaN	94.240900	16.579494	0.732306	0.181167	0.378867	2023-01
6	2023-01-01 01:00:00	M1	A	298.584877	14.081397	78.070456	32.066631	81.757595	78.382944	102.908197	...	NaN	NaN	NaN	NaN	122.067285	16.796991	0.768684	1.529528	0.985470	2023-01
7	2023-01-01 01:00:00	M2	A	263.573060	13.701868	75.966263	30.772568	68.424088	89.181910	103.415857	...	NaN	NaN	NaN	NaN	101.313846	16.857293	0.722970	0.390418	0.154186	2023-01
8	2023-01-01 01:00:00	M3	A	285.552436	13.476033	74.714179	29.745000	69.694977	84.289027	103.312593	...	NaN	NaN	NaN	NaN	103.550331	19.181352	0.660425	0.315991	0.116196	2023-01
9	2023-01-01 01:00:00	M4	A	306.647052	13.862049	76.854344	31.962362	71.580560	73.980311	108.120992	...	NaN	NaN	NaN	NaN	110.339358	16.717480	0.768506	0.180361	0.050174	2023-01

10 rows × 65 columns

2. Análisis Univariado

2.1 Variables Categóricas

# Análisis de variables categóricas
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
fig.suptitle('Distribución de Variables Categóricas', fontsize=16, fontweight='bold')

# Molino ID
molino_counts = df['molino_id'].value_counts()
axes[0,0].bar(molino_counts.index, molino_counts.values, color='steelblue', alpha=0.7)
axes[0,0].set_title('Distribución por Molino')
axes[0,0].set_xlabel('Molino ID')
axes[0,0].set_ylabel('Frecuencia')

# Turno
turno_counts = df['turno'].value_counts()
axes[0,1].bar(turno_counts.index, turno_counts.values, color='orange', alpha=0.7)
axes[0,1].set_title('Distribución por Turno')
axes[0,1].set_xlabel('Turno')
axes[0,1].set_ylabel('Frecuencia')

# Tipo de falla
falla_counts = df['tipo_falla'].value_counts()
axes[1,0].bar(range(len(falla_counts)), falla_counts.values, color='red', alpha=0.7)
axes[1,0].set_title('Distribución por Tipo de Falla')
axes[1,0].set_xlabel('Tipo de Falla')
axes[1,0].set_ylabel('Frecuencia')
axes[1,0].set_xticks(range(len(falla_counts)))
axes[1,0].set_xticklabels(falla_counts.index, rotation=45)

# Severidad de falla
severidad_counts = df['severidad_falla'].value_counts().sort_index()
axes[1,1].bar(severidad_counts.index, severidad_counts.values, color='purple', alpha=0.7)
axes[1,1].set_title('Distribución por Severidad de Falla')
axes[1,1].set_xlabel('Severidad')
axes[1,1].set_ylabel('Frecuencia')

plt.tight_layout()
plt.show()

# Estadísticas de variables categóricas
print("=== ESTADÍSTICAS VARIABLES CATEGÓRICAS ===")
print(f"Molinos únicos: {df['molino_id'].nunique()}")
print(f"Turnos únicos: {df['turno'].nunique()}")
print(f"Tipos de falla únicos: {df['tipo_falla'].nunique()}")
print(f"Niveles de severidad: {df['severidad_falla'].nunique()}")

=== ESTADÍSTICAS VARIABLES CATEGÓRICAS ===
Molinos únicos: 6
Turnos únicos: 3
Tipos de falla únicos: 5
Niveles de severidad: 4

2.2 Variables de Proceso y Operación

# Definición de grupos de variables
variables_proceso = [
    'feed_rate', 'velocidad_rotacion', 'velocidad_porcentaje_critica',
    'nivel_carga_bolas', 'densidad_pulpa', 'agua_adicionada', 'presion_ciclones'
]

variables_monitoreo = [
    'vibracion_cojinete_feed_h', 'vibracion_cojinete_feed_v',
    'vibracion_cojinete_discharge_h', 'vibracion_cojinete_discharge_v',
    'temp_cojinete_feed', 'temp_cojinete_discharge', 'corriente_motor'
]

variables_resultado = [
    'throughput_real', 'eficiencia_molienda', 'consumo_energetico_especifico',
    'granulometria_producto_p80', 'eficiencia_clasificacion'
]

variables_mineral = [
    'work_index_bond', 'dureza_mineral', 'humedad_mineral',
    'granulometria_feed_p80', 'densidad_mineral', 'contenido_arcillas'
]

# Función para crear histogramas con estadísticas
def plot_histograms_with_stats(variables, data, title, rows=2, cols=4):
    n_vars = len(variables)
    fig, axes = plt.subplots(rows, cols, figsize=(20, 10))
    fig.suptitle(title, fontsize=16, fontweight='bold')
    
    axes = axes.flatten() if rows * cols > 1 else [axes]
    
    for i, var in enumerate(variables):
        if i < len(axes):
            ax = axes[i]
            values = data[var].dropna()
            
            # Histograma
            ax.hist(values, bins=30, alpha=0.7, color='steelblue', edgecolor='black')
            ax.set_title(f'{var}')
            ax.set_xlabel('Valor')
            ax.set_ylabel('Frecuencia')
            
            # Estadísticas en el gráfico
            mean_val = values.mean()
            std_val = values.std()
            ax.axvline(mean_val, color='red', linestyle='--', alpha=0.8, label=f'Media: {mean_val:.2f}')
            ax.axvline(mean_val + std_val, color='orange', linestyle=':', alpha=0.8)
            ax.axvline(mean_val - std_val, color='orange', linestyle=':', alpha=0.8)
            ax.legend(fontsize=8)
            
            # Grid
            ax.grid(True, alpha=0.3)
    
    # Ocultar ejes sobrantes
    for i in range(len(variables), len(axes)):
        axes[i].set_visible(False)
    
    plt.tight_layout()
    plt.show()
    
    # Tabla de estadísticas
    stats_df = data[variables].describe().round(2)
    return stats_df

# Análisis de variables de proceso
print("=== VARIABLES DE PROCESO ===")
stats_proceso = plot_histograms_with_stats(variables_proceso, df, 'Distribución Variables de Proceso')
print(stats_proceso)

=== VARIABLES DE PROCESO ===

       feed_rate  velocidad_rotacion  velocidad_porcentaje_critica  \
count  131478.00           131478.00                     131478.00   
mean      278.69               13.71                         76.01   
std        23.97                0.36                          1.99   
min       154.71               12.63                         70.00   
25%       262.46               13.47                         74.67   
50%       278.35               13.71                         76.00   
75%       294.62               13.95                         77.35   
max       406.25               15.13                         83.87   

       nivel_carga_bolas  densidad_pulpa  agua_adicionada  presion_ciclones  
count          131478.00       131478.00        131478.00         131478.00  
mean               32.00           72.07            86.66             95.07  
std                 1.49            4.26            13.79             15.49  
min                28.00           29.94            42.88             28.46  
25%                30.99           69.36            76.67             84.49  
50%                32.01           72.04            86.26             94.92  
75%                33.02           74.80            96.36            105.46  
max                36.00           90.92           138.87            164.05  

2.3 Detección de Outliers

# Función para detectar outliers usando IQR
def detect_outliers_iqr(data, columns):
    outliers_summary = {}
    
    for col in columns:
        Q1 = data[col].quantile(0.25)
        Q3 = data[col].quantile(0.75)
        IQR = Q3 - Q1
        lower_bound = Q1 - 1.5 * IQR
        upper_bound = Q3 + 1.5 * IQR
        
        outliers = data[(data[col] < lower_bound) | (data[col] > upper_bound)]
        outliers_summary[col] = {
            'count': len(outliers),
            'percentage': (len(outliers) / len(data)) * 100,
            'lower_bound': lower_bound,
            'upper_bound': upper_bound
        }
    
    return outliers_summary

# Detectar outliers en variables numéricas clave
numeric_columns = df.select_dtypes(include=[np.number]).columns.tolist()
numeric_columns = [col for col in numeric_columns if col not in ['horas_operacion_acumuladas', 'ciclos_arranque_parada', 'severidad_falla']]

outliers_info = detect_outliers_iqr(df, numeric_columns[:15])  # Primeras 15 variables numéricas

# Crear DataFrame con información de outliers
outliers_df = pd.DataFrame(outliers_info).T
outliers_df = outliers_df.sort_values('percentage', ascending=False)

print("=== DETECCIÓN DE OUTLIERS (Top 10) ===")
print(outliers_df.head(10))

# Boxplots para variables con más outliers
top_outlier_vars = outliers_df.head(6).index.tolist()

fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Boxplots - Variables con Mayor Presencia de Outliers', fontsize=16, fontweight='bold')

axes = axes.flatten()
for i, var in enumerate(top_outlier_vars):
    df.boxplot(column=var, ax=axes[i])
    axes[i].set_title(f'{var}\n({outliers_df.loc[var, "percentage"]:.1f}% outliers)')
    axes[i].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

=== DETECCIÓN DE OUTLIERS (Top 10) ===
                                  count  percentage  lower_bound  upper_bound
vibracion_cojinete_discharge_h  11433.0    8.695751     4.479105    11.938371
vibracion_cojinete_feed_v       11410.0    8.678258     3.513618    10.150157
vibracion_cojinete_discharge_v  11294.0    8.590030     4.034109    11.569205
vibracion_cojinete_feed_h       11106.0    8.447041     3.902710    10.437621
vibracion_pinion                 1200.0    0.912700     4.930961    10.440259
densidad_pulpa                   1153.0    0.876953    61.196987    82.968277
feed_rate                        1059.0    0.805458   214.213435   342.864201
vibracion_shell_v                 968.0    0.736245     4.394161     8.444470
vibracion_shell_h                 928.0    0.705822     4.384596     8.434432
velocidad_rotacion                904.0    0.687568    12.743850    14.677068

3.Análisis de Fallas y Mantenimiento Predictivo

3.1 Distribución Temporal de Fallas

# Análisis de fallas por tiempo
print("=== ANÁLISIS DE FALLAS ===")

# Distribución de fallas en diferentes horizontes
fallas_summary = {
    'Fallas en 7 días': df['falla_en_7d'].sum(),
    'Fallas en 14 días': df['falla_en_14d'].sum(),
    'Fallas en 30 días': df['falla_en_30d'].sum()
}

print("Distribución de fallas por horizonte temporal:")
for key, value in fallas_summary.items():
    percentage = (value / len(df)) * 100
    print(f"  {key}: {value} casos ({percentage:.2f}%)")

# Visualización de fallas por molino y tipo
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('Análisis de Fallas por Molino y Tipo', fontsize=16, fontweight='bold')

# Fallas por molino
fallas_molino = df.groupby('molino_id')[['falla_en_7d', 'falla_en_14d', 'falla_en_30d']].sum()
fallas_molino.plot(kind='bar', ax=axes[0,0], color=['red', 'orange', 'yellow'])
axes[0,0].set_title('Fallas por Molino')
axes[0,0].set_xlabel('Molino ID')
axes[0,0].set_ylabel('Número de Fallas')
axes[0,0].legend(['7 días', '14 días', '30 días'])
axes[0,0].tick_params(axis='x', rotation=0)

# Distribución de tipos de falla
tipo_falla_dist = df['tipo_falla'].value_counts()
axes[0,1].pie(tipo_falla_dist.values, labels=tipo_falla_dist.index, autopct='%1.1f%%')
axes[0,1].set_title('Distribución por Tipo de Falla')

# Severidad de fallas por molino
severidad_molino = df.groupby(['molino_id', 'severidad_falla']).size().unstack(fill_value=0)
severidad_molino.plot(kind='bar', stacked=True, ax=axes[1,0], colormap='Reds')
axes[1,0].set_title('Severidad de Fallas por Molino')
axes[1,0].set_xlabel('Molino ID')
axes[1,0].set_ylabel('Número de Casos')
axes[1,0].tick_params(axis='x', rotation=0)

# Distribución de días hasta falla
dias_falla = df[df['dias_hasta_falla'] < 365]['dias_hasta_falla']
axes[1,1].hist(dias_falla, bins=30, alpha=0.7, color='red', edgecolor='black')
axes[1,1].set_title('Distribución de Días hasta Falla')
axes[1,1].set_xlabel('Días hasta Falla')
axes[1,1].set_ylabel('Frecuencia')
axes[1,1].axvline(dias_falla.mean(), color='blue', linestyle='--', label=f'Media: {dias_falla.mean():.1f}')
axes[1,1].legend()

plt.tight_layout()
plt.show()

=== ANÁLISIS DE FALLAS ===
Distribución de fallas por horizonte temporal:
  Fallas en 7 días: 5544 casos (4.22%)
  Fallas en 14 días: 11088 casos (8.43%)
  Fallas en 30 días: 23412 casos (17.81%)

3.2 Análisis de Señales Precursoras

# Análisis de variables que pueden predecir fallas
def analyze_failure_predictors(data, failure_col, variables_to_analyze):
    """Analiza variables que pueden predecir fallas"""
    
    results = {}
    
    for var in variables_to_analyze:
        # Separar datos con y sin falla
        no_failure = data[data[failure_col] == False][var]
        with_failure = data[data[failure_col] == True][var]
        
        # Test estadístico
        statistic, p_value = stats.mannwhitneyu(no_failure.dropna(), with_failure.dropna(), alternative='two-sided')
        
        # Diferencia de medias
        mean_diff = with_failure.mean() - no_failure.mean()
        mean_diff_pct = (mean_diff / no_failure.mean()) * 100 if no_failure.mean() != 0 else 0
        
        results[var] = {
            'mean_no_failure': no_failure.mean(),
            'mean_with_failure': with_failure.mean(),
            'mean_difference': mean_diff,
            'mean_diff_percentage': mean_diff_pct,
            'p_value': p_value,
            'significant': p_value < 0.05
        }
    
    return pd.DataFrame(results).T

# Variables a analizar como predictores
predictor_vars = [
    'vibracion_cojinete_feed_h', 'vibracion_cojinete_feed_v',
    'temp_cojinete_feed', 'temp_cojinete_discharge', 'corriente_motor',
    'consumo_energetico_especifico', 'eficiencia_molienda',
    'nivel_desgaste_liners', 'calidad_aceite_ppm',
    'anomaly_score_vibration', 'anomaly_score_electrical'
]

# Análisis para fallas en 7 días
print("=== ANÁLISIS DE PREDICTORES DE FALLAS (7 días) ===")
predictors_7d = analyze_failure_predictors(df, 'falla_en_7d', predictor_vars)
predictors_7d_sig = predictors_7d[predictors_7d['significant']].sort_values('mean_diff_percentage', key=abs, ascending=False)

print("Variables significativamente diferentes entre casos con y sin falla:")
print(predictors_7d_sig[['mean_no_failure', 'mean_with_failure', 'mean_diff_percentage', 'p_value']].round(4))

# Visualización de predictores más significativos
top_predictors = predictors_7d_sig.head(6).index.tolist()

fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Distribución de Variables Predictoras (Fallas en 7 días)', fontsize=16, fontweight='bold')

axes = axes.flatten()
for i, var in enumerate(top_predictors):
    # Datos para boxplot
    no_failure_data = df[df['falla_en_7d'] == False][var]
    with_failure_data = df[df['falla_en_7d'] == True][var]
    
    # Boxplot comparativo
    data_to_plot = [no_failure_data.dropna(), with_failure_data.dropna()]
    box_plot = axes[i].boxplot(data_to_plot, labels=['Sin Falla', 'Con Falla'], patch_artist=True)
    
    # Colorear boxplots
    box_plot['boxes'][0].set_facecolor('lightblue')
    box_plot['boxes'][1].set_facecolor('lightcoral')
    
    axes[i].set_title(f'{var}')
    axes[i].grid(True, alpha=0.3)
    
    # Añadir información estadística
    diff_pct = predictors_7d_sig.loc[var, 'mean_diff_percentage']
    p_val = predictors_7d_sig.loc[var, 'p_value']
    axes[i].text(0.5, 0.95, f'Diff: {diff_pct:.1f}%\np-val: {p_val:.3f}', 
                transform=axes[i].transAxes, ha='center', va='top',
                bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

plt.tight_layout()
plt.show()

=== ANÁLISIS DE PREDICTORES DE FALLAS (7 días) ===
Variables significativamente diferentes entre casos con y sin falla:
                          mean_no_failure mean_with_failure  \
anomaly_score_vibration           0.36147          3.625876   
vibracion_cojinete_feed_v        7.099313         22.289988   
vibracion_cojinete_feed_h        7.436884         22.788696   
anomaly_score_electrical         0.805707          0.841328   
temp_cojinete_feed              70.430238          72.12789   
eficiencia_molienda             74.725947          74.54171   

                          mean_diff_percentage   p_value  
anomaly_score_vibration             903.090762       0.0  
vibracion_cojinete_feed_v           213.973884       0.0  
vibracion_cojinete_feed_h           206.428029       0.0  
anomaly_score_electrical              4.421152  0.000003  
temp_cojinete_feed                    2.410402       0.0  
eficiencia_molienda                  -0.246551  0.000682  

4. Análisis Bivariado y Correlaciones

4.1 Matriz de Correlación por Grupos

# Matriz de correlación para variables de proceso
def plot_correlation_matrix(data, variables, title, figsize=(12, 10)):
    """Crea matriz de correlación con anotaciones"""
    corr_matrix = data[variables].corr()
    
    plt.figure(figsize=figsize)
    mask = np.triu(np.ones_like(corr_matrix, dtype=bool))
    
    sns.heatmap(corr_matrix, mask=mask, annot=True, cmap='RdBu_r', center=0,
                square=True, linewidths=0.5, cbar_kws={"shrink": 0.8}, fmt='.2f')
    
    plt.title(title, fontsize=14, fontweight='bold')
    plt.tight_layout()
    plt.show()
    
    return corr_matrix

# Correlaciones por grupos de variables
print("=== MATRICES DE CORRELACIÓN POR GRUPOS ===")

# Variables de proceso
corr_proceso = plot_correlation_matrix(df, variables_proceso, 'Correlaciones - Variables de Proceso')

# Variables de monitoreo
corr_monitoreo = plot_correlation_matrix(df, variables_monitoreo, 'Correlaciones - Variables de Monitoreo')

# Variables de resultado
corr_resultado = plot_correlation_matrix(df, variables_resultado, 'Correlaciones - Variables de Resultado')

=== MATRICES DE CORRELACIÓN POR GRUPOS ===

4.2 Relaciones Clave: Eficiencia vs Variables Operacionales

# Análisis de relaciones con eficiencia
print("=== ANÁLISIS DE EFICIENCIA ===")

# Variables clave que afectan la eficiencia
efficiency_vars = [
    'feed_rate', 'velocidad_porcentaje_critica', 'nivel_carga_bolas',
    'densidad_pulpa', 'presion_ciclones', 'work_index_bond'
]

# Correlaciones con eficiencia de molienda
efficiency_corr = df[efficiency_vars + ['eficiencia_molienda']].corr()['eficiencia_molienda'].drop('eficiencia_molienda')
efficiency_corr_sorted = efficiency_corr.abs().sort_values(ascending=False)

print("Correlaciones con Eficiencia de Molienda:")
for var in efficiency_corr_sorted.index:
    corr_val = efficiency_corr[var]
    print(f"  {var}: {corr_val:.3f}")

# Scatter plots de variables más correlacionadas
top_efficiency_vars = efficiency_corr_sorted.head(4).index.tolist()

fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('Relaciones con Eficiencia de Molienda', fontsize=16, fontweight='bold')

axes = axes.flatten()
for i, var in enumerate(top_efficiency_vars):
    scatter = axes[i].scatter(df[var], df['eficiencia_molienda'], 
                            c=df['molino_id'].astype('category').cat.codes, 
                            alpha=0.6, cmap='tab10')
    axes[i].set_xlabel(var)
    axes[i].set_ylabel('Eficiencia Molienda (%)')
    axes[i].set_title(f'Eficiencia vs {var}')
    axes[i].grid(True, alpha=0.3)
    
    # Línea de tendencia
    z = np.polyfit(df[var].dropna(), df.loc[df[var].notna(), 'eficiencia_molienda'], 1)
    p = np.poly1d(z)
    axes[i].plot(df[var], p(df[var]), "r--", alpha=0.8, linewidth=2)
    
    # Correlación en el gráfico
    corr_val = efficiency_corr[var]
    axes[i].text(0.05, 0.95, f'r = {corr_val:.3f}', transform=axes[i].transAxes,
                bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

plt.tight_layout()
plt.show()

=== ANÁLISIS DE EFICIENCIA ===
Correlaciones con Eficiencia de Molienda:
  densidad_pulpa: 0.007
  presion_ciclones: 0.003
  feed_rate: -0.003
  nivel_carga_bolas: -0.003
  velocidad_porcentaje_critica: 0.002
  work_index_bond: -0.000

4.3 Consumo Energético vs Variables de Control

# Análisis del consumo energético
print("\n=== ANÁLISIS DE CONSUMO ENERGÉTICO ===")

energy_vars = [
    'feed_rate', 'velocidad_rotacion', 'nivel_carga_bolas',
    'densidad_pulpa', 'dureza_mineral', 'work_index_bond', 'throughput_real'
]

# Correlaciones con consumo energético específico
energy_corr = df[energy_vars + ['consumo_energetico_especifico']].corr()['consumo_energetico_especifico'].drop('consumo_energetico_especifico')
energy_corr_sorted = energy_corr.abs().sort_values(ascending=False)

print("Correlaciones con Consumo Energético Específico:")
for var in energy_corr_sorted.index:
    corr_val = energy_corr[var]
    print(f"  {var}: {corr_val:.3f}")

# Análisis por molino
energy_by_molino = df.groupby('molino_id')['consumo_energetico_especifico'].agg(['mean', 'std', 'count'])
print("\nConsumo energético por molino:")
print(energy_by_molino.round(2))

# Visualización del consumo energético
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('Análisis de Consumo Energético', fontsize=16, fontweight='bold')

# Distribución por molino
df.boxplot(column='consumo_energetico_especifico', by='molino_id', ax=axes[0,0])
axes[0,0].set_title('Consumo Energético por Molino')
axes[0,0].set_xlabel('Molino ID')
axes[0,0].set_ylabel('Consumo Energético Específico (kWh/t)')

# Consumo vs Throughput
scatter = axes[0,1].scatter(df['throughput_real'], df['consumo_energetico_especifico'],
                          c=df['eficiencia_molienda'], cmap='viridis', alpha=0.6)
axes[0,1].set_xlabel('Throughput Real (t/h)')
axes[0,1].set_ylabel('Consumo Energético Específico (kWh/t)')
axes[0,1].set_title('Consumo vs Throughput (Color: Eficiencia)')
plt.colorbar(scatter, ax=axes[0,1], label='Eficiencia Molienda (%)')

# Consumo vs Work Index
axes[1,0].scatter(df['work_index_bond'], df['consumo_energetico_especifico'],
                 alpha=0.6, color='orange')
axes[1,0].set_xlabel('Work Index Bond (kWh/t)')
axes[1,0].set_ylabel('Consumo Energético Específico (kWh/t)')
axes[1,0].set_title('Consumo vs Work Index')
axes[1,0].grid(True, alpha=0.3)

# Tendencia temporal del consumo
df_monthly = df.groupby(df['timestamp'].dt.to_period('M'))['consumo_energetico_especifico'].mean()
axes[1,1].plot(df_monthly.index.astype(str), df_monthly.values, 'bo-', linewidth=2, markersize=6)
axes[1,1].set_xlabel('Período')
axes[1,1].set_ylabel('Consumo Promedio (kWh/t)')
axes[1,1].set_title('Tendencia Temporal del Consumo')
axes[1,1].tick_params(axis='x', rotation=45)
axes[1,1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

=== ANÁLISIS DE CONSUMO ENERGÉTICO ===
Correlaciones con Consumo Energético Específico:
  work_index_bond: 0.828
  dureza_mineral: 0.506
  densidad_pulpa: 0.083
  feed_rate: 0.062
  nivel_carga_bolas: 0.003
  velocidad_rotacion: 0.003
  throughput_real: 0.000

Consumo energético por molino:
            mean   std  count
molino_id                    
M1         15.90  1.44  21913
M2         15.29  1.40  21913
M3         16.59  1.51  21913
M4         15.44  1.40  21913
M5         16.23  1.47  21913
M6         15.75  1.43  21913

5. Análisis Multivariado

5.1 Análisis de Componentes Principales (PCA)

# Preparación de datos para PCA
print("=== ANÁLISIS DE COMPONENTES PRINCIPALES ===")

# Seleccionar variables numéricas para PCA (excluyendo variables target y timestamps)
pca_vars = [col for col in df.select_dtypes(include=[np.number]).columns 
           if col not in ['horas_operacion_acumuladas', 'ciclos_arranque_parada', 
                         'severidad_falla', 'dias_hasta_falla']]

# Preparar datos
X_pca = df[pca_vars].dropna()
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_pca)

# Aplicar PCA
pca = PCA()
X_pca_transformed = pca.fit_transform(X_scaled)

# Análisis de varianza explicada
cumsum_var = np.cumsum(pca.explained_variance_ratio_)
n_components_90 = np.argmax(cumsum_var >= 0.90) + 1
n_components_95 = np.argmax(cumsum_var >= 0.95) + 1

print(f"Componentes para explicar 90% de varianza: {n_components_90}")
print(f"Componentes para explicar 95% de varianza: {n_components_95}")

# Visualización del PCA
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('Análisis de Componentes Principales', fontsize=16, fontweight='bold')

# Varianza explicada
axes[0,0].plot(range(1, min(21, len(pca.explained_variance_ratio_)+1)), 
               pca.explained_variance_ratio_[:20], 'bo-', linewidth=2, markersize=6)
axes[0,0].set_xlabel('Componente Principal')
axes[0,0].set_ylabel('Varianza Explicada')
axes[0,0].set_title('Varianza Explicada por Componente')
axes[0,0].grid(True, alpha=0.3)

# Varianza acumulada
axes[0,1].plot(range(1, min(21, len(cumsum_var)+1)), cumsum_var[:20], 'ro-', linewidth=2, markersize=6)
axes[0,1].axhline(y=0.90, color='green', linestyle='--', alpha=0.8, label='90%')
axes[0,1].axhline(y=0.95, color='orange', linestyle='--', alpha=0.8, label='95%')
axes[0,1].set_xlabel('Número de Componentes')
axes[0,1].set_ylabel('Varianza Explicada Acumulada')
axes[0,1].set_title('Varianza Acumulada')
axes[0,1].legend()
axes[0,1].grid(True, alpha=0.3)

# Scatter plot PC1 vs PC2 por molino
molino_colors = df.loc[X_pca.index, 'molino_id']
for i, molino in enumerate(df['molino_id'].unique()):
    mask = molino_colors == molino
    axes[1,0].scatter(X_pca_transformed[mask, 0], X_pca_transformed[mask, 1], 
                     alpha=0.6, label=f'Molino {molino}', s=30)

axes[1,0].set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.1%} varianza)')
axes[1,0].set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.1%} varianza)')
axes[1,0].set_title('Proyección en PC1 vs PC2 por Molino')
axes[1,0].legend()
axes[1,0].grid(True, alpha=0.3)

# Loadings de las primeras componentes
feature_names = pca_vars
pc1_loadings = pca.components_[0]
pc2_loadings = pca.components_[1]

# Top features por componente
top_pc1_idx = np.argsort(np.abs(pc1_loadings))[-8:]
top_pc2_idx = np.argsort(np.abs(pc2_loadings))[-8:]

axes[1,1].barh(range(len(top_pc1_idx)), pc1_loadings[top_pc1_idx], alpha=0.7, color='blue', label='PC1')
axes[1,1].barh(range(len(top_pc2_idx)), pc2_loadings[top_pc2_idx], alpha=0.7, color='red', label='PC2')
axes[1,1].set_yticks(range(max(len(top_pc1_idx), len(top_pc2_idx))))
axes[1,1].set_yticklabels([feature_names[i] for i in top_pc1_idx], fontsize=8)
axes[1,1].set_xlabel('Loading')
axes[1,1].set_title('Top Loadings PC1 y PC2')
axes[1,1].legend()

plt.tight_layout()
plt.show()

# Top loadings para PC1 y PC2
print("\nTop Variables en PC1:")
pc1_contributions = list(zip(feature_names, pc1_loadings))
pc1_sorted = sorted(pc1_contributions, key=lambda x: abs(x[1]), reverse=True)
for var, loading in pc1_sorted[:10]:
    print(f"  {var}: {loading:.3f}")

print("\nTop Variables en PC2:")
pc2_contributions = list(zip(feature_names, pc2_loadings))
pc2_sorted = sorted(pc2_contributions, key=lambda x: abs(x[1]), reverse=True)
for var, loading in pc2_sorted[:10]:
    print(f"  {var}: {loading:.3f}")

=== ANÁLISIS DE COMPONENTES PRINCIPALES ===
Componentes para explicar 90% de varianza: 27
Componentes para explicar 95% de varianza: 31

Top Variables en PC1:
  temp_cojinete_feed: 0.275
  temp_cojinete_discharge: 0.275
  potencia_activa: 0.265
  corriente_motor: 0.265
  consumo_energetico_especifico: 0.258
  temp_motor_principal: 0.257
  work_index_bond: 0.250
  potencia_especifica_neta: 0.250
  temperatura_trend_7d: 0.248
  energia_trend_24h: 0.240

Top Variables en PC2:
  anomaly_score_vibration: 0.419
  vibracion_cojinete_discharge_v: 0.402
  vibracion_cojinete_feed_h: 0.401
  vibracion_cojinete_discharge_h: 0.401
  vibracion_cojinete_feed_v: 0.401
  vibracion_trend_7d: 0.330
  consumo_energetico_especifico: -0.095
  potencia_especifica_neta: -0.092
  energia_trend_24h: -0.092
  eficiencia_molienda: 0.081

6. Análisis Temporal

6.1 Tendencias Temporales

# Análisis de tendencias temporales
print("=== ANÁLISIS TEMPORAL ===")

# Preparar datos temporales
df['date'] = df['timestamp'].dt.date
df['hour'] = df['timestamp'].dt.hour
df['day_of_week'] = df['timestamp'].dt.day_name()
df['month'] = df['timestamp'].dt.month_name()

# Tendencias mensuales de variables clave
monthly_trends = df.groupby(df['timestamp'].dt.to_period('M')).agg({
    'eficiencia_molienda': 'mean',
    'consumo_energetico_especifico': 'mean',
    'throughput_real': 'mean',
    'vibracion_cojinete_feed_h': 'mean',
    'temp_cojinete_feed': 'mean',
    'falla_en_7d': 'sum'
}).round(2)

# Visualización de tendencias
fig, axes = plt.subplots(3, 2, figsize=(16, 18))
fig.suptitle('Análisis de Tendencias Temporales', fontsize=16, fontweight='bold')

# Eficiencia mensual
axes[0,0].plot(monthly_trends.index.astype(str), monthly_trends['eficiencia_molienda'], 
               'go-', linewidth=2, markersize=6)
axes[0,0].set_title('Tendencia - Eficiencia de Molienda')
axes[0,0].set_ylabel('Eficiencia (%)')
axes[0,0].tick_params(axis='x', rotation=45)
axes[0,0].grid(True, alpha=0.3)

# Consumo energético mensual
axes[0,1].plot(monthly_trends.index.astype(str), monthly_trends['consumo_energetico_especifico'], 
               'ro-', linewidth=2, markersize=6)
axes[0,1].set_title('Tendencia - Consumo Energético')
axes[0,1].set_ylabel('Consumo (kWh/t)')
axes[0,1].tick_params(axis='x', rotation=45)
axes[0,1].grid(True, alpha=0.3)

# Throughput mensual
axes[1,0].plot(monthly_trends.index.astype(str), monthly_trends['throughput_real'], 
               'bo-', linewidth=2, markersize=6)
axes[1,0].set_title('Tendencia - Throughput Real')
axes[1,0].set_ylabel('Throughput (t/h)')
axes[1,0].tick_params(axis='x', rotation=45)
axes[1,0].grid(True, alpha=0.3)

# Vibración mensual
axes[1,1].plot(monthly_trends.index.astype(str), monthly_trends['vibracion_cojinete_feed_h'], 
               'mo-', linewidth=2, markersize=6)
axes[1,1].set_title('Tendencia - Vibración Cojinete')
axes[1,1].set_ylabel('Vibración (mm/s)')
axes[1,1].tick_params(axis='x', rotation=45)
axes[1,1].grid(True, alpha=0.3)

# Temperatura mensual
axes[2,0].plot(monthly_trends.index.astype(str), monthly_trends['temp_cojinete_feed'], 
               'co-', linewidth=2, markersize=6)
axes[2,0].set_title('Tendencia - Temperatura Cojinete')
axes[2,0].set_ylabel('Temperatura (°C)')
axes[2,0].tick_params(axis='x', rotation=45)
axes[2,0].grid(True, alpha=0.3)

# Fallas mensuales
axes[2,1].bar(monthly_trends.index.astype(str), monthly_trends['falla_en_7d'], 
              alpha=0.7, color='red')
axes[2,1].set_title('Fallas por Mes')
axes[2,1].set_ylabel('Número de Fallas')
axes[2,1].tick_params(axis='x', rotation=45)
axes[2,1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("Tendencias mensuales:")
print(monthly_trends)

=== ANÁLISIS TEMPORAL ===

Tendencias mensuales:
           eficiencia_molienda  consumo_energetico_especifico  \
timestamp                                                       
2023-01                  74.73                          16.08   
2023-02                  74.69                          16.80   
2023-03                  74.73                          17.17   
2023-04                  74.73                          17.22   
2023-05                  74.77                          16.78   
2023-06                  74.70                          16.11   
2023-07                  74.67                          15.27   
2023-08                  74.72                          14.53   
2023-09                  74.65                          14.14   
2023-10                  74.75                          14.09   
2023-11                  74.72                          14.52   
2023-12                  74.79                          15.22   
2024-01                  74.72                          16.08   
2024-02                  74.71                          16.79   
2024-03                  74.69                          17.21   
2024-04                  74.75                          17.21   
2024-05                  74.76                          16.81   
2024-06                  74.73                          16.07   
2024-07                  74.78                          15.28   
2024-08                  74.74                          14.53   
2024-09                  74.74                          14.09   
2024-10                  74.71                          14.11   
2024-11                  74.72                          14.55   
2024-12                  74.68                          15.28   
2025-01                  74.71                          16.11   
2025-02                  74.67                          16.79   
2025-03                  74.70                          17.19   
2025-04                  74.73                          17.21   
2025-05                  74.63                          16.84   
2025-06                  74.74                          16.11   
2025-07                  74.36                          15.77   

           throughput_real  vibracion_cojinete_feed_h  temp_cojinete_feed  \
timestamp                                                                   
2023-01             264.09                       8.16               71.28   
2023-02             264.82                       7.61               74.24   
2023-03             264.69                       9.71               75.99   
2023-04             264.67                       8.01               75.96   
2023-05             264.07                      11.33               74.51   
2023-06             264.57                       8.10               71.59   
2023-07             264.80                       6.66               67.90   
2023-08             264.50                       6.50               64.93   
2023-09             264.38                       7.04               63.36   
2023-10             263.94                       8.22               63.41   
2023-11             264.91                       7.16               64.99   
2023-12             264.54                       6.65               67.84   
2024-01             264.43                       7.13               71.14   
2024-02             264.66                      10.60               74.55   
2024-03             264.81                       8.69               75.91   
2024-04             264.85                       9.14               76.06   
2024-05             264.26                       7.92               74.17   
2024-06             264.67                       8.32               71.38   
2024-07             264.71                       7.42               68.02   
2024-08             264.43                       8.26               65.19   
2024-09             264.56                       7.19               63.30   
2024-10             264.08                       8.84               63.61   
2024-11             264.21                       6.60               65.05   
2024-12             264.60                       7.36               68.09   
2025-01             264.73                       7.13               71.33   
2025-02             264.56                       8.66               74.35   
2025-03             264.27                       8.44               76.20   
2025-04             264.51                      11.29               76.16   
2025-05             264.06                       7.47               74.17   
2025-06             263.91                       7.08               71.19   
2025-07             263.95                       6.99               70.17   

           falla_en_7d  
timestamp               
2023-01            168  
2023-02              0  
2023-03            336  
2023-04              0  
2023-05            503  
2023-06            337  
2023-07              0  
2023-08              0  
2023-09            168  
2023-10            336  
2023-11            168  
2023-12              0  
2024-01              0  
2024-02            336  
2024-03            168  
2024-04            168  
2024-05             42  
2024-06            153  
2024-07            141  
2024-08            504  
2024-09            168  
2024-10            504  
2024-11              0  
2024-12            336  
2025-01              0  
2025-02            168  
2025-03            336  
2025-04            504  
2025-05              0  
2025-06              0  
2025-07              0  

6.2 Análisis por Turnos

# Análisis por turnos
print("\n=== ANÁLISIS POR TURNOS ===")

# Estadísticas por turno
turno_stats = df.groupby('turno').agg({
    'eficiencia_molienda': ['mean', 'std'],
    'consumo_energetico_especifico': ['mean', 'std'],
    'throughput_real': ['mean', 'std'],
    'falla_en_7d': 'sum',
    'temp_cojinete_feed': 'mean',
    'vibracion_cojinete_feed_h': 'mean'
}).round(2)

print("Estadísticas por turno:")
print(turno_stats)

# Visualización por turnos
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Análisis por Turnos de Trabajo', fontsize=16, fontweight='bold')

# Eficiencia por turno
df.boxplot(column='eficiencia_molienda', by='turno', ax=axes[0,0])
axes[0,0].set_title('Eficiencia por Turno')
axes[0,0].set_ylabel('Eficiencia (%)')

# Consumo por turno
df.boxplot(column='consumo_energetico_especifico', by='turno', ax=axes[0,1])
axes[0,1].set_title('Consumo Energético por Turno')
axes[0,1].set_ylabel('Consumo (kWh/t)')

# Throughput por turno
df.boxplot(column='throughput_real', by='turno', ax=axes[0,2])
axes[0,2].set_title('Throughput por Turno')
axes[0,2].set_ylabel('Throughput (t/h)')

# Temperatura por turno
df.boxplot(column='temp_cojinete_feed', by='turno', ax=axes[1,0])
axes[1,0].set_title('Temperatura por Turno')
axes[1,0].set_ylabel('Temperatura (°C)')

# Vibración por turno
df.boxplot(column='vibracion_cojinete_feed_h', by='turno', ax=axes[1,1])
axes[1,1].set_title('Vibración por Turno')
axes[1,1].set_ylabel('Vibración (mm/s)')

# Fallas por turno
fallas_turno = df.groupby('turno')['falla_en_7d'].sum()
axes[1,2].bar(fallas_turno.index, fallas_turno.values, alpha=0.7, color='red')
axes[1,2].set_title('Fallas por Turno')
axes[1,2].set_ylabel('Número de Fallas')

plt.tight_layout()
plt.show()

# Test estadístico entre turnos
from scipy.stats import kruskal

print("\nTests estadísticos entre turnos (Kruskal-Wallis):")
variables_test = ['eficiencia_molienda', 'consumo_energetico_especifico', 'throughput_real']

for var in variables_test:
    groups = [df[df['turno'] == turno][var].dropna() for turno in df['turno'].unique()]
    statistic, p_value = kruskal(*groups)
    significance = "Significativo" if p_value < 0.05 else "No significativo"
    print(f"  {var}: p-value = {p_value:.4f} ({significance})")

=== ANÁLISIS POR TURNOS ===
Estadísticas por turno:
      eficiencia_molienda       consumo_energetico_especifico        \
                     mean   std                          mean   std   
turno                                                                 
A                   74.72  3.94                         15.86  1.51   
B                   74.72  3.95                         15.87  1.51   
C                   74.71  3.96                         15.87  1.51   

      throughput_real        falla_en_7d temp_cojinete_feed  \
                 mean    std         sum               mean   
turno                                                         
A              258.20  20.19        2079              69.66   
B              270.06  19.95        1848              71.23   
C              266.16  19.84        1617              70.75   

      vibracion_cojinete_feed_h  
                           mean  
turno                            
A                          7.89  
B                          8.26  
C                          8.13  

Tests estadísticos entre turnos (Kruskal-Wallis):
  eficiencia_molienda: p-value = 0.9366 (No significativo)
  consumo_energetico_especifico: p-value = 0.5754 (No significativo)
  throughput_real: p-value = 0.0000 (Significativo)

6.3 Patrones Horarios y Estacionales

# Análisis de patrones horarios
print("\n=== PATRONES HORARIOS ===")

# Promedio por hora del día
hourly_patterns = df.groupby('hour').agg({
    'eficiencia_molienda': 'mean',
    'consumo_energetico_especifico': 'mean',
    'temp_cojinete_feed': 'mean',
    'throughput_real': 'mean'
}).round(2)

# Visualización de patrones horarios
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('Patrones Horarios de Operación', fontsize=16, fontweight='bold')

# Eficiencia por hora
axes[0,0].plot(hourly_patterns.index, hourly_patterns['eficiencia_molienda'], 
               'go-', linewidth=2, markersize=6)
axes[0,0].set_title('Eficiencia por Hora del Día')
axes[0,0].set_xlabel('Hora')
axes[0,0].set_ylabel('Eficiencia (%)')
axes[0,0].grid(True, alpha=0.3)
axes[0,0].set_xticks(range(0, 24, 2))

# Consumo por hora
axes[0,1].plot(hourly_patterns.index, hourly_patterns['consumo_energetico_especifico'], 
               'ro-', linewidth=2, markersize=6)
axes[0,1].set_title('Consumo Energético por Hora')
axes[0,1].set_xlabel('Hora')
axes[0,1].set_ylabel('Consumo (kWh/t)')
axes[0,1].grid(True, alpha=0.3)
axes[0,1].set_xticks(range(0, 24, 2))

# Temperatura por hora
axes[1,0].plot(hourly_patterns.index, hourly_patterns['temp_cojinete_feed'], 
               'co-', linewidth=2, markersize=6)
axes[1,0].set_title('Temperatura por Hora del Día')
axes[1,0].set_xlabel('Hora')
axes[1,0].set_ylabel('Temperatura (°C)')
axes[1,0].grid(True, alpha=0.3)
axes[1,0].set_xticks(range(0, 24, 2))

# Throughput por hora
axes[1,1].plot(hourly_patterns.index, hourly_patterns['throughput_real'], 
               'bo-', linewidth=2, markersize=6)
axes[1,1].set_title('Throughput por Hora del Día')
axes[1,1].set_xlabel('Hora')
axes[1,1].set_ylabel('Throughput (t/h)')
axes[1,1].grid(True, alpha=0.3)
axes[1,1].set_xticks(range(0, 24, 2))

plt.tight_layout()
plt.show()

print("Patrones horarios promedio:")
print(hourly_patterns)

=== PATRONES HORARIOS ===

Patrones horarios promedio:
      eficiencia_molienda  consumo_energetico_especifico  temp_cojinete_feed  \
hour                                                                           
0                   74.71                          15.84               69.45   
1                   74.77                          15.82               69.43   
2                   74.79                          15.85               69.37   
3                   74.78                          15.85               69.48   
4                   74.71                          15.89               69.55   
5                   74.70                          15.86               69.44   
6                   74.62                          15.91               69.53   
7                   74.72                          15.88               69.46   
8                   74.66                          15.85               71.27   
9                   74.68                          15.87               71.38   
10                  74.72                          15.86               71.29   
11                  74.72                          15.88               71.31   
12                  74.72                          15.88               71.33   
13                  74.74                          15.86               71.27   
14                  74.75                          15.86               71.27   
15                  74.73                          15.88               71.36   
16                  74.73                          15.84               70.61   
17                  74.79                          15.86               70.66   
18                  74.60                          15.90               70.78   
19                  74.74                          15.86               70.73   
20                  74.72                          15.87               70.83   
21                  74.81                          15.84               70.61   
22                  74.65                          15.91               70.81   
23                  74.70                          15.89               70.82   

      throughput_real  
hour                   
0              256.72  
1              257.14  
2              256.41  
3              256.82  
4              256.66  
5              256.67  
6              256.30  
7              256.47  
8              270.57  
9              270.76  
10             270.68  
11             270.21  
12             270.83  
13             270.53  
14             270.88  
15             270.79  
16             265.81  
17             265.96  
18             266.10  
19             266.43  
20             266.38  
21             265.82  
22             266.04  
23             266.38  

7. Análisis de Anomalías

# Análisis de anomalías usando scores calculados
print("=== ANÁLISIS DE ANOMALÍAS ===")

# Distribución de scores de anomalía
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('Análisis de Scores de Anomalía', fontsize=16, fontweight='bold')

# Distribución de anomaly scores
axes[0,0].hist(df['anomaly_score_vibration'], bins=30, alpha=0.7, color='red', edgecolor='black')
axes[0,0].set_title('Distribución - Anomaly Score Vibración')
axes[0,0].set_xlabel('Score')
axes[0,0].set_ylabel('Frecuencia')
axes[0,0].axvline(df['anomaly_score_vibration'].mean(), color='blue', linestyle='--', 
                 label=f'Media: {df["anomaly_score_vibration"].mean():.3f}')
axes[0,0].legend()

axes[0,1].hist(df['anomaly_score_electrical'], bins=30, alpha=0.7, color='orange', edgecolor='black')
axes[0,1].set_title('Distribución - Anomaly Score Eléctrico')
axes[0,1].set_xlabel('Score')
axes[0,1].set_ylabel('Frecuencia')
axes[0,1].axvline(df['anomaly_score_electrical'].mean(), color='blue', linestyle='--',
                 label=f'Media: {df["anomaly_score_electrical"].mean():.3f}')
axes[0,1].legend()

# Relación entre anomaly scores y fallas
no_falla = df[df['falla_en_7d'] == False]
con_falla = df[df['falla_en_7d'] == True]

axes[1,0].scatter(no_falla['anomaly_score_vibration'], no_falla['anomaly_score_electrical'], 
                 alpha=0.6, color='blue', label='Sin Falla', s=20)
axes[1,0].scatter(con_falla['anomaly_score_vibration'], con_falla['anomaly_score_electrical'], 
                 alpha=0.8, color='red', label='Con Falla', s=30, marker='^')
axes[1,0].set_xlabel('Anomaly Score Vibración')
axes[1,0].set_ylabel('Anomaly Score Eléctrico')
axes[1,0].set_title('Relación Anomaly Scores vs Fallas')
axes[1,0].legend()
axes[1,0].grid(True, alpha=0.3)

# Evolución temporal de anomaly scores
monthly_anomalies = df.groupby(df['timestamp'].dt.to_period('M')).agg({
    'anomaly_score_vibration': 'mean',
    'anomaly_score_electrical': 'mean'
})

axes[1,1].plot(monthly_anomalies.index.astype(str), monthly_anomalies['anomaly_score_vibration'], 
               'ro-', linewidth=2, markersize=6, label='Vibración')
axes[1,1].plot(monthly_anomalies.index.astype(str), monthly_anomalies['anomaly_score_electrical'], 
               'bo-', linewidth=2, markersize=6, label='Eléctrico')
axes[1,1].set_title('Evolución Temporal Anomaly Scores')
axes[1,1].set_ylabel('Score Promedio')
axes[1,1].legend()
axes[1,1].tick_params(axis='x', rotation=45)
axes[1,1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Estadísticas de anomalías por molino
print("Anomaly scores promedio por molino:")
anomaly_by_molino = df.groupby('molino_id')[['anomaly_score_vibration', 'anomaly_score_electrical']].mean()
print(anomaly_by_molino.round(4))

# Correlación entre anomaly scores y variables de monitoreo
anomaly_correlations = df[['anomaly_score_vibration', 'anomaly_score_electrical', 
                          'vibracion_cojinete_feed_h', 'temp_cojinete_feed', 
                          'corriente_motor', 'eficiencia_molienda']].corr()

print("\nCorrelaciones con variables de monitoreo:")
print(anomaly_correlations[['anomaly_score_vibration', 'anomaly_score_electrical']].round(3))

=== ANÁLISIS DE ANOMALÍAS ===

Anomaly scores promedio por molino:
           anomaly_score_vibration  anomaly_score_electrical
molino_id                                                   
M1                          0.5385                    0.8057
M2                          0.5155                    0.8066
M3                          0.4491                    0.8086
M4                          0.5143                    0.8076
M5                          0.4822                    0.8072
M6                          0.4951                    0.8075

Correlaciones con variables de monitoreo:
                           anomaly_score_vibration  anomaly_score_electrical
anomaly_score_vibration                      1.000                     0.071
anomaly_score_electrical                     0.071                     1.000
vibracion_cojinete_feed_h                    0.875                     0.039
temp_cojinete_feed                          -0.029                     0.033
corriente_motor                             -0.089                     0.077
eficiencia_molienda                          0.040                    -0.044

8. Insights y Recomendaciones

8.1 Hallazgos Principales

1. Análisis de eficiencia por molino

eficiencia_molino = df.groupby('molino_id')['eficiencia_molienda'].agg(['mean', 'std'])
mejor_molino = eficiencia_molino['mean'].idxmax()
peor_molino = eficiencia_molino['mean'].idxmin()
print(f"\n1. RENDIMIENTO POR MOLINO:")
print(f"   • Mejor rendimiento: {mejor_molino} ({eficiencia_molino.loc[mejor_molino, 'mean']:.1f}% eficiencia)")
print(f"   • Menor rendimiento: {peor_molino} ({eficiencia_molino.loc[peor_molino, 'mean']:.1f}% eficiencia)")
print(f"   • Diferencia: {eficiencia_molino.loc[mejor_molino, 'mean'] - eficiencia_molino.loc[peor_molino, 'mean']:.1f} puntos porcentuales")

1. RENDIMIENTO POR MOLINO:
   • Mejor rendimiento: M2 (77.5% eficiencia)
   • Menor rendimiento: M3 (71.4% eficiencia)
   • Diferencia: 6.1 puntos porcentuales

2. Análisis de consumo energético

consumo_molino = df.groupby('molino_id')['consumo_energetico_especifico'].mean()
molino_eficiente = consumo_molino.idxmin()
molino_ineficiente = consumo_molino.idxmax()
print(f"\n2. CONSUMO ENERGÉTICO:")
print(f"   • Más eficiente: {molino_eficiente} ({consumo_molino[molino_eficiente]:.1f} kWh/t)")
print(f"   • Menos eficiente: {molino_ineficiente} ({consumo_molino[molino_ineficiente]:.1f} kWh/t)")
print(f"   • Diferencia: {consumo_molino[molino_ineficiente] - consumo_molino[molino_eficiente]:.1f} kWh/t")

2. CONSUMO ENERGÉTICO:
   • Más eficiente: M2 (15.3 kWh/t)
   • Menos eficiente: M3 (16.6 kWh/t)
   • Diferencia: 1.3 kWh/t

3. Análisis de fallas

fallas_molino = df.groupby('molino_id')['falla_en_7d'].sum()
molino_problematico = fallas_molino.idxmax()
molino_confiable = fallas_molino.idxmin()
print(f"\n3. CONFIABILIDAD:")
print(f"   • Más problemático: {molino_problematico} ({fallas_molino[molino_problematico]} fallas)")
print(f"   • Más confiable: {molino_confiable} ({fallas_molino[molino_confiable]} fallas)")

3. CONFIABILIDAD:
   • Más problemático: M1 (1008 fallas)
   • Más confiable: M2 (840 fallas)

4. Variables críticas para predicción

if 'predictors_7d_sig' in locals():
    top_predictors_list = predictors_7d_sig.head(5).index.tolist()
    print(f"\n4. VARIABLES PREDICTORAS CRÍTICAS:")
    for i, var in enumerate(top_predictors_list, 1):
        diff_pct = predictors_7d_sig.loc[var, 'mean_diff_percentage']
        print(f"   {i}. {var}: {diff_pct:+.1f}% diferencia pre-falla")

4. VARIABLES PREDICTORAS CRÍTICAS:
   1. anomaly_score_vibration: +903.1% diferencia pre-falla
   2. vibracion_cojinete_feed_v: +214.0% diferencia pre-falla
   3. vibracion_cojinete_feed_h: +206.4% diferencia pre-falla
   4. anomaly_score_electrical: +4.4% diferencia pre-falla
   5. temp_cojinete_feed: +2.4% diferencia pre-falla

5. Optimización energética

energy_efficiency_corr = df[['work_index_bond', 'dureza_mineral', 'throughput_real', 'consumo_energetico_especifico']].corr()['consumo_energetico_especifico']
print(f"\n5. FACTORES DE CONSUMO ENERGÉTICO:")
print(f"   • Work Index Bond: r = {energy_efficiency_corr['work_index_bond']:.3f}")
print(f"   • Dureza Mineral: r = {energy_efficiency_corr['dureza_mineral']:.3f}")
print(f"   • Throughput: r = {energy_efficiency_corr['throughput_real']:.3f}")

5. FACTORES DE CONSUMO ENERGÉTICO:
   • Work Index Bond: r = 0.828
   • Dureza Mineral: r = 0.506
   • Throughput: r = 0.000

6. Patrones temporales

turno_mejor = df.groupby('turno')['eficiencia_molienda'].mean().idxmax()
turno_peor = df.groupby('turno')['eficiencia_molienda'].mean().idxmin()
print(f"\n6. PATRONES OPERACIONALES:")
print(f"   • Mejor turno: {turno_mejor}")
print(f"   • Turno con mayor consumo: {df.groupby('turno')['consumo_energetico_especifico'].mean().idxmax()}")

6. PATRONES OPERACIONALES:
   • Mejor turno: B
   • Turno con mayor consumo: C

8.2 Dashboard de KPIs Críticos

# Dashboard de métricas clave
fig = make_subplots(
    rows=3, cols=3,
    subplot_titles=[
        'Eficiencia por Molino', 'Consumo Energético por Molino', 'Fallas por Molino',
        'Distribución de Eficiencia', 'Tendencia Mensual Eficiencia', 'Correlación Eficiencia-Consumo',
        'Anomaly Scores', 'Disponibilidad por Turno', 'ROI Potencial'
    ],
    specs=[[{"type": "bar"}, {"type": "bar"}, {"type": "bar"}],
           [{"type": "histogram"}, {"type": "scatter"}, {"type": "scatter"}],
           [{"type": "scatter"}, {"type": "bar"}, {"type": "bar"}]]
)

# Row 1: Métricas por molino
molinos = df['molino_id'].unique()
eficiencia_avg = [df[df['molino_id'] == m]['eficiencia_molienda'].mean() for m in molinos]
consumo_avg = [df[df['molino_id'] == m]['consumo_energetico_especifico'].mean() for m in molinos]
fallas_count = [df[df['molino_id'] == m]['falla_en_7d'].sum() for m in molinos]

fig.add_trace(go.Bar(x=molinos, y=eficiencia_avg, name='Eficiencia', 
                     marker_color='green'), row=1, col=1)
fig.add_trace(go.Bar(x=molinos, y=consumo_avg, name='Consumo', 
                     marker_color='red'), row=1, col=2)
fig.add_trace(go.Bar(x=molinos, y=fallas_count, name='Fallas', 
                     marker_color='orange'), row=1, col=3)

# Row 2: Distribuciones y tendencias
fig.add_trace(go.Histogram(x=df['eficiencia_molienda'], name='Dist. Eficiencia',
                          marker_color='lightblue'), row=2, col=1)

# Tendencia mensual
monthly_eff = df.groupby(df['timestamp'].dt.to_period('M'))['eficiencia_molienda'].mean()
fig.add_trace(go.Scatter(x=monthly_eff.index.astype(str), y=monthly_eff.values,
                        mode='lines+markers', name='Tendencia', 
                        line=dict(color='blue')), row=2, col=2)

# Correlación eficiencia-consumo
fig.add_trace(go.Scatter(x=df['eficiencia_molienda'], y=df['consumo_energetico_especifico'],
                        mode='markers', name='Eff vs Consumo',
                        marker=dict(color='purple', size=4, opacity=0.6)), row=2, col=3)

# Row 3: Anomalías y análisis operacional
fig.add_trace(go.Scatter(x=df['anomaly_score_vibration'], y=df['anomaly_score_electrical'],
                        mode='markers', name='Anomalías',
                        marker=dict(color='red', size=4, opacity=0.6)), row=3, col=1)

# Disponibilidad por turno (eficiencia como proxy)
turno_stats = df.groupby('turno')['eficiencia_molienda'].mean()
fig.add_trace(go.Bar(x=turno_stats.index, y=turno_stats.values, 
                     name='Eficiencia por Turno', marker_color='cyan'), row=3, col=2)

# ROI potencial (basado en diferencias de eficiencia)
roi_data = []
for molino in molinos:
    current_eff = df[df['molino_id'] == molino]['eficiencia_molienda'].mean()
    potential_improvement = max(eficiencia_avg) - current_eff
    roi_data.append(potential_improvement)

fig.add_trace(go.Bar(x=molinos, y=roi_data, name='Mejora Potencial (%)', 
                     marker_color='gold'), row=3, col=3)

# Actualizar layout
fig.update_layout(
    height=1200,
    title_text="Dashboard de KPIs Críticos - Molinos MineraPeru",
    title_x=0.5,
    showlegend=False,
    font=dict(size=10)
)

fig.show()

8.3 Matriz de Priorización de Acciones

# Matriz de priorización basada en impacto y facilidad de implementación
print("\n=== MATRIZ DE PRIORIZACIÓN DE ACCIONES ===")

acciones = {
    'Acción': [
        'Optimizar parámetros Molino M6',
        'Mantenimiento predictivo vibraciones',
        'Ajuste alimentación por Work Index',
        'Estandarizar mejores prácticas Turno A',
        'Monitoreo temperatura cojinetes',
        'Optimizar densidad pulpa',
        'Calibración sensores Molino M5',
        'Reducir variabilidad operacional',
        'Implementar control automático',
        'Training operadores clusters alto rendimiento'
    ],
    'Impacto_Financiero': [9, 8, 7, 6, 8, 7, 5, 8, 9, 6],  # 1-10
    'Facilidad_Implementacion': [8, 6, 7, 9, 7, 8, 9, 5, 3, 8],  # 1-10
    'Molino_Objetivo': ['M6', 'Todos', 'Todos', 'Todos', 'M5,M6', 'Todos', 'M5', 'Todos', 'Todos', 'Todos'],
    'Tiempo_Implementacion': ['2 sem', '4 sem', '3 sem', '1 sem', '2 sem', '2 sem', '1 sem', '8 sem', '12 sem', '4 sem']
}

matriz_acciones = pd.DataFrame(acciones)
matriz_acciones['Score_Prioridad'] = matriz_acciones['Impacto_Financiero'] * matriz_acciones['Facilidad_Implementacion']
matriz_acciones = matriz_acciones.sort_values('Score_Prioridad', ascending=False)

print("Acciones priorizadas por impacto y facilidad:")
print(matriz_acciones[['Acción', 'Score_Prioridad', 'Molino_Objetivo', 'Tiempo_Implementacion']])

# Visualización de la matriz
fig, ax = plt.subplots(figsize=(12, 8))
scatter = ax.scatter(matriz_acciones['Facilidad_Implementacion'], 
                    matriz_acciones['Impacto_Financiero'],
                    s=matriz_acciones['Score_Prioridad']*3,
                    alpha=0.6, c=range(len(matriz_acciones)), cmap='viridis')

# Añadir etiquetas
for i, txt in enumerate(matriz_acciones.index):
    ax.annotate(f'{txt+1}', (matriz_acciones.iloc[i]['Facilidad_Implementacion'], 
                            matriz_acciones.iloc[i]['Impacto_Financiero']),
                fontsize=10, ha='center', va='center')

ax.set_xlabel('Facilidad de Implementación (1-10)')
ax.set_ylabel('Impacto Financiero (1-10)')
ax.set_title('Matriz de Priorización de Acciones\n(Tamaño = Score de Prioridad)')
ax.grid(True, alpha=0.3)

# Líneas de cuadrantes
ax.axhline(y=6.5, color='red', linestyle='--', alpha=0.5)
ax.axvline(x=6.5, color='red', linestyle='--', alpha=0.5)

# Etiquetas de cuadrantes
ax.text(8.5, 8.5, 'GANAR RÁPIDO\n(Alta prioridad)', ha='center', va='center', 
        bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.7))
ax.text(3.5, 8.5, 'PROYECTOS GRANDES\n(Planificar)', ha='center', va='center',
        bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.7))
ax.text(8.5, 3.5, 'LLENAR ESPACIOS\n(Si hay tiempo)', ha='center', va='center',
        bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.7))
ax.text(3.5, 3.5, 'TRABAJO DURO\n(Baja prioridad)', ha='center', va='center',
        bbox=dict(boxstyle='round', facecolor='lightcoral', alpha=0.7))

plt.tight_layout()
plt.show()

=== MATRIZ DE PRIORIZACIÓN DE ACCIONES ===
Acciones priorizadas por impacto y facilidad:
                                          Acción  Score_Prioridad  \
0                 Optimizar parámetros Molino M6               72   
4                Monitoreo temperatura cojinetes               56   
5                       Optimizar densidad pulpa               56   
3         Estandarizar mejores prácticas Turno A               54   
2             Ajuste alimentación por Work Index               49   
1           Mantenimiento predictivo vibraciones               48   
9  Training operadores clusters alto rendimiento               48   
6                 Calibración sensores Molino M5               45   
7               Reducir variabilidad operacional               40   
8                 Implementar control automático               27   

  Molino_Objetivo Tiempo_Implementacion  
0              M6                 2 sem  
4           M5,M6                 2 sem  
5           Todos                 2 sem  
3           Todos                 1 sem  
2           Todos                 3 sem  
1           Todos                 4 sem  
9           Todos                 4 sem  
6              M5                 1 sem  
7           Todos                 8 sem  
8           Todos                12 sem  

9. Conclusiones y Recomendaciones Ejecutivas

9.1 Conclusiones Principales

🎯 Oportunidades de Mejora Identificadas:

1. Disparidad entre Molinos: Existe una diferencia significativa en eficiencia entre molinos, con el mejor superando al peor en varios puntos porcentuales, representando una oportunidad inmediata de mejora.
2. Predictores de Fallas: Se identificaron variables clave que muestran comportamiento anómalo días antes de las fallas, especialmente en vibraciones y temperaturas.
3. Optimización Energética: El consumo energético específico muestra alta correlación con características del mineral y parámetros operacionales, sugiriendo oportunidades de optimización dinámica.
4. Patrones Operacionales: Existen diferencias sistemáticas entre turnos y patrones horarios que indican oportunidades de estandarización de mejores prácticas.

9.2 Recomendaciones Estratégicas

🚀 Acciones Inmediatas (0-3 meses):

• Implementar sistema de alertas tempranas basado en variables predictoras identificadas
• Estandarizar parámetros operacionales del molino de mejor rendimiento
• Optimizar schedules de mantenimiento basado en análisis de supervivencia

📈 Iniciativas Mediano Plazo (3-12 meses):

• Desarrollar modelo de optimización dinámica para ajuste de parámetros según características del mineral
• Implementar sistema de monitoreo continuo de anomalías
• Establecer programa de mejora continua basado en clustering operacional

🔬 Proyectos Largo Plazo (12+ meses):

• Automatización de control de proceso basado en machine learning
• Integración completa de sistema predictivo con mantenimiento
• Expansión del modelo a otros circuitos de la planta

9.3 Impacto Financiero Estimado

Basado en los hallazgos del análisis, las mejoras propuestas podrían generar:

• Reducción de costos operativos: 15-20% mediante optimización energética
• Incremento en throughput: 8-12% por mayor disponibilidad
• Reducción costos mantenimiento: 25-30% por mantenimiento predictivo
• ROI estimado: 300-400% en el primer año

💡 Próximos Pasos:

1. Validar hallazgos con equipo de operaciones
2. Desarrollar prototipo de sistema predictivo
3. Diseñar programa piloto en molino seleccionado
4. Establecer métricas de seguimiento y KPIs de éxito